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cn ! Abstract 

I propose a method to calculate logarithmic interaction in two dimensions and coulomb inter- 
action in three dimensions under periodic boundary conditions. This paper considers the case of 



a rectangular cell in two dimensions and an orthorhombic cell in three dimensions. Unlike the 
Ewald method 1 , there is no parameter to be optimized, nor does it involve error functions, thus 
leading to the accuracy obtained. This method is similar in approach to that of Sperb^ , but the 
derivation is considerably simpler and physically appealing. An important aspect of the proposed 
method is the faster convergence of the Green function for a particular case as compared to Sperb's 
work. The convergence of the sums for the most part of unit cell is exponential, and hence requires 
the calculation of only a few dozen terms. In a very simple way, we also obtain expressions for 

> 

t^j- . interaction for systems with slab geometries. Expressions for the Madelung constant of CsCl and 

o . 

I/"") , NaCl are also obtained. 
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I. INTRODUCTION 



In Molecular dynamics (MD) and Monte Carlo (MC) simulations one is required to cal- 
culate the potential energy and forces acting on a particle due to other particles. Sometimes 
such forces have a long range interaction. In such situations, periodic boundary condi- 
tions are usually imposed in order to avoid the boundary effects, which might be especially 
prominent for small systems that are usually employed in MD simulations. Under periodic 
boundary conditions interaction of a particle with another particle includes the direct inter- 
action plus an interaction of the first particle with all replicas of itself as well as all replicas 
of the second particle. These replicas come into picture due to the periodic repetitions of a 
charge under the periodic boundary conditions. The energy contribution arising from the 
interaction of a particle with its own replicas is termed as the self energy. The calculation 
of self energy is important in an MC simulation, where size of the box might change during 
simulation, such as in isobaric MC. The natural question that arises is how may one compute 
the long range interaction of a particle with a second particle along with all the replicas of 
the second particle. The self energy part may then be obtained trivially as well. For eighty 
years, researchers have employed the Ewald sum technique to perform such summations. 
However, the Ewald sum technique has certain drawbacks. The primary drawback being, 
the optimization of a parameter that renders break up of the original algebraic sum in two 
parts, one in real space and the other one in Fourier space. Only when this parameter is 
chosen properly do the sums in real and Fourier spaces converge fast. A second problem 
with the Ewald sum is that even if one achieves optimal choice of the parameter for break- 
ing up the sum, one might lose numerical accuracy as the terms in these sums involve error 
functions, whose evaluation to high degree of accuracy is difficult. In this paper we will 
consider the logarithmic interaction in two dimensions and Coulomb interaction in three 
di—, The 2D ease has bee, satisfactorily dea lt with in Re, Q Thus m ainl y we will 
concentrate on 3D results. The Ewald method is the most widely used technique for system 
in 3D. An alternative technique for summation over long range forces in 3D for a cubic 
unit cell was given by Lekner—. A tedious method was employed to obtained the self energy 
part of the interaction. However, Lekner generalized his work to an orthorhombic cell^ and 
obtained self energies in a much simpler manner. These recent methods by Lekner- and 
Sperb^ are similar in spirit but their derivation involves complicated algebra. One problem 
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with Lekner's expressions is that they involve a triple sum. Sperb's 2 - results are better in 
that part of the interaction has only a double sum. Nevertheless a triple sum (Eq. 2.4 and 
2.7 in Ref. 0) is still employed for the case when both particles are very close to each other. 

The technique that we propose is based on a series summation in Fourier space. Work 
along these lines has been previously reported in recent papers^, as well as by Harris et 
al&, SperVj 2 -, Crandall et a/.— and Marshall^. The outline of this paper follows. In section 
II, we derive a general formula for dimension d > 2. In section III the formula is applied 
to get logarithmic sum in 2D. Section IV describes application of the general formula to 
get Coulomb summation 1/r for the slab geometry well as for 3D case. Section 

V considers evaluation of Madelung constants for CsCl and NaCl. Finally, we discuss our 
results in Section VI. 



II. COULOMB SUM IN d DIMENSION 



An interaction in which satisfies the Poisson equation in d dimensions will be termed 
as a Coulomb type potential for that particular dimension. For example the logarithmic 
interaction is a Coulomb type interaction in 2D. In this section we discuss how one can 
calculate a pairwise Coulomb interaction between two particles, separated by a displacement 
r. For simplicity, we consider the case of a unit charge situated within an orthorhombic cell 
in d dimensions. Let the d sides of the unit cell be labeled by Zi, The basic unit cell 

repeats itself in all d dimensions. The unit charge interacts with other identical unit charges 
(for the case of different charges qi and qi one just gets an extra factor of q\q-i) situated at 
the vertices of the periodic structure. The interaction between two particles is given by the 
Green's function in d dimension, G(r), which satisfies the Poisson equation, 

V 2 G(r) = -C d ^5(r + l). (2.1) 
l 

where V 2 is the Laplacian operator in d dimensions, 1 denotes a d dimensional vector, whose 
components are integer multiples of Zj's and Cd is specified by 

_ [ B 2 for d = 2, 

I (d - 2) B d for d > 2, 
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where B d stands for the coefficient of {d — 1) dimensional surface element in d dimensions, 



d in) 2 

B d = f ) ! - . 2.2 

r(f+i) v ; 

Here T(x) stands for the gamma function. Thus B 2 = 2tt, B% = 4n etc. We note that coef- 
ficients in Eq. (|2.1|) have been chosen such that G(r) stands for a Coulomb type summation 
in d dimensions. For example, for d > 2, we will have G given by, 

G(x h x 2 ,--- ,Xd) = ^2 ~ O"' ( 2 - 3 ) 

{m} d {(™>ih ~ xi) 2 + (m 2 l 2 ~ x 2 ) 2 H h (m d l d - x d ) 2 } 2 

where {m} d stands for set of d numbers mi,m 2 , ...,m d . The summation over each runs 
from — oo to +oo. The solution to Eq. (j2.1j) can be easily expressed in Fourier space, 

G( Xl , x 2 ,-.., x d ) = r, 7-2 . 7-2 T^Y' ( 2 - 4 ) 

[ } 12 d ^{{ir) + 

where < Xi/k < 1. The function G(xi,x 2 ,- ■■ ,x d ), as defined above, diverges since the 
term corresponding to all m's being equal to zero blows up. This is expected since the sum 
defined in Eq. (l2.4|) has contribution coming from an infinite set of identical charges. For the 
sum in Eq. f!2.4|) to make sense we add an infinitesimal term to the denominator and subtract 
off a counter term from the whole sum as follows: 



G(x ly x 2 , ••■ ,x d ) 



a, 
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(27T) 2 /!/ 2 
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(2.5) 
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J2TT(m 1 ^-+m 2 ^--\ hm d |^) 
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where £ is an infinitesimal parameter which tends to zero. The prescription employed above 
amounts to assumption of the presence of a uniform background charge. For example, let 
us consider the case of 3D. For every charge, q, one may imagine a uniform distribution of 
charge, such that, total charge per unit cell adds up to —q. For a charge neutral periodic 
system, imposing these kind of background uniform charge distributions does not matter 
since total uniform background charge adds up to zero. However, now a unit charge located 
within the unit cell at position (xi,x 2 ,xs) not only interacts with a second charge located 



at the origin and its periodic images, but it also interacts with the neutralizing background 
charge, compensating the charge of the second particle. This particular way of introducing 
the neutralizing background charge leads to only the intrinsic part4 of potential energy. 
Now, it can be easily verified that Eq. (|2.5|) satisfies the following equation: 

V 2 G(r) = -C d S(r + 1) + (2.6) 

where the last term in Eq. (J2.6|) represents the uniform background charge. Complete ex- 
pression for the potential has a term arising from surface contribution. For the 2D case this 
turns out to be zero, but for 3D one obtains a contribution from dipole ternA 
Moving further, we can perform one of the d sums in Eq.flgfl analytically^-, 

oo i2nm d ^ 

9(Xd ' {mU) = w + kw 2 +•••+' + e (27) 



n cosh 



7r 7d ({m},e) (l-2Jal 



ld({m},£) sinh[7T7 d ({m},0] 
where stands for li/lj and 7d({ m },0 is defined as 



ld({m} , = yj(m 1 l d>1 ) 2 + --- + (m d _ 1 l d>d _ 1 ) 2 + e. (2.8) 
For convenience we also define 

ldo{{m} , = ^ (miZ d ,i) 2 + • • • + (m d _iZ did _i) 2 . (2.9) 
Using Eqs. (|2.5j) and ()2.7|) one obtains 

G{x 1 ,x 2 ,--- ,x d ) = . \ ld — (2.10) 



x 



( tt ( xA 1 

\{m}d-i i=1 1 



In the limit £ — > 0, the term corresponding to all being set to zero in Eq. (|2.1()|) must be 
separated out as follows: 
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G(xi,x 2 , ■ ■ ■ 




(2-11) 



where a prime on the summation sign implies that the term corresponding to all m; being 
zero is not to be included. In Eq. ([2.11jh we separated out the term corresponding to all m ; 
being set to zero and took the limit £ — > as follows, 



/ T cosh k(l- 2^)1 j\ ^ ( /, ,x / yi 

n [i U M - e J = t I 1 - 6 (V) + 6 U) } • (2l2) 

Eq. ()2.1H) forms the main result derived in this section. It is important to note that as a 
result of the symmetry present in the problem, it suffices to look at only that part of the 
unit cell which corresponds to < Xi/U < 0.5 for all i's. Hence, from here on we will assume 
< Xi/U < 0.5. In the next two sections, we investigate two important cases corresponding 
to d = 2 and d = 3. 



III. LOGARITHMIC SUM IN TWO DIMENSIONS 



Energy of N particles contained in a rectangular unit cell with periodic boundaries and 
interacting through a logarithmic potential in 2D can be expressed as^, 

^totai = \ E WG*1?< - r i) + E 

where charges are denoted by qi and the position of charges in the unit cell by where 
1 < i < N. We will obtain expressions for G?2d(r) and G^ lf in this section. The pairwise 
interaction is given by the Green function G2d( r ) which satisfies the Poisson equation in 2D, 

2-7T 

V 2 G 2d (r) = -27T V 5(r + 1) + — , (3.2) 
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where the last term on the rhs of Eq. (j3.2j) stands for the neutralizing background charge. 
Eq. (|3.2j) is a special case of Eq. (j2.1jl . We look for a solution of Eq. (j3.2|) with periodic 
boundary conditions along X\ and £2 directions. This solution can be easily obtained from 
the general formula, Eq. (J2.11)) . derived in the previous section, 



G2d(%l, %2 




2 M 



sinh [7T72 (m)] 
2' 

I .To \ 

+ 6 



— COS I 2,71171 — 

h 



(3.3) 




where a prime on m implies the term corresponding to m = is to be excluded. Without 
any loss of generality we may assume that sides of the rhombic cells have been labeled so 
that l\ < 12- This condition will make sure that 720 ( m ) > 1 f° r a H integer values of m. Let 
us now consider the convergence of the sum in Eq. (J3.3J) . The first part of Eq. (j3.3|) converges 
exponentially, but in some cases the convergence may be very slow. Specifically, the leading 
term in ()3.3|) decays as exp (— 27r|m| \x2\/h))- Thus the convergence depends on the ratio 
x%jl\. We see that one obtains a slow exponential convergence when < x%jl\ < 0.1. To 
handle this case properly, we break the first sum in Eq. (j3.3|) into two parts by application 
of a trigonometric identity, 



cosh(a — b) cosh(a) exp(— b) 



sinh (6) 

This leads to the expression: 



2tt ^ 



1 \ ^ 7T 



sinh (b) 



exp (—a J 



7T 


cosh 


7rm/2,i 
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\m 


sinh(7r 


m 


hi) 



cos lum 



+ 



-E 

7r L — ' m 

m=l 
1 X 

TT 5^ 



exp (— 7r \m\ /2,i) cosh 



irm l 2 ,i ( 2f 



7T 



7T * — ' m 

m=l 



exp ( —2nm 



sinh(7rm / 2j i 
'x 2 



COS ( ZTTm — 

h 



cos 2nm 



(3.4) 



(3.5) 



We notice that the first part of Eq. (J3.5|) converges even for the case when < X2/I1 < 0.1. 
In fact the slowest convergence for the first part will now occur for the case when 2x2 = h- 
But even this "slowest" convergence amounts to a very rapid exponential convergence of 
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exp (— 7i\m\l 2 /li)). We have yet to account for the last sum in Eq. (j3.5J) . Using the identity 

OO j 

— exp(— 2nnx) cos(2irny) = (3-6) 



— In < cosh f 27r^ | — cos f 2tt^- 



7T- 



,x > 0, 



n=l 

1 In (2) 
— — In [cosh(27ra;) — cos(27r?/)J + nx — 

the last part of the sum in Eq. (|3.5jl may be explicitly evaluated to 

] f - (3.7) 

Assembling the terms together, we finally obtain the following expression for the 2D Green 
function, 



2tt 



\ In {cosh [2^ 
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sinh 
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cos 2nm 



cos I 2ir— 
V h 

In (2) 



Self energy may be easily obtained as 



^2d 
^self 



lim t G 2d (xi,x 2 ) + In I Jxj + x\ 

(xi,a;2)^(0,0) ' 



1 



E 



i: exp^-7r|m|| 



-In 



2nZ ^f M sinh(7r|m||^ V h 



2ty 



7T l 2 

6 h' 



(3.8) 



(3.9) 



The results derived here may be trivially generalized to the case of a rhombic cell, but our 
concern in this paper has only been with orthorhombic cases. The results obtained here 
were numerically checked and found to be in agreement with those of Gr0nbech- Jensen^.. 

IV. COULOMB SUM IN 3D 



Energy of N particles contained in a orthorhombic unit cell with periodic boundaries and 
interacting through a Coulomb type potential in 3D can be expressed as, 



3d , 2tt i 

self ^ q I 2-^i ™ 



(4.1) 



where charges are denoted by g« and the position of charges in the unit cell by and 
1 < i < N . We will obtain expressions for G3d(r) and G^if 111 this section. The application 
of Eq. (J2.11j) for an orthorhombic cell in 3D leads to 



1 U 71 



cosh 



*73o(M) (l - 2^) 



G 3d (x l ,x 2 ,x 3 ) = --j- V n —- —r-7 77 U1 (4.2) 



X 



n( x i\ ^3 TT j 




where 



7so({m} , £) = yjirmhrf + (m 2 / 3 , 2 ) 2 . (4.3) 
Without any loss of generality we assume that axis have been labeled such that 

h > h > h- (4.4) 

The condition in Eq. (j4.4j) makes sure 730 ({ m }) > 1 for all sets {m}. Eq. (|4.2|) is 
one of our main results for 3D case. We note that the potential energy obtained 
consists of only the intrinsic part 4 . A dipole contribution will have to be included 
in Eq. (J4.2|) to obtain the real potential energy 4 "^. This dipole contribution is repre- 
sented by the last term on the rhs in Eq. (j4.1|) We notice that the sum in Eq. fl4.2j) con- 
verges exponentially. In fact the terms corresponding to large \mi\ and \m,2\ decay as 
exp ^—2nx 3 ^ (mi/Ii) 2 + {m 2 /l 2 ) 2 ^j , which with the assumption in Eq. (j4.4j) means that 

terms decay faster than exp (^—2txx 3 \J (mi/l 2 ) 2 + (rn 2 /l 2 ) 2 ^j ■ Thus the convergence depends 
upon the ratio r 32 = x 3 /l 2 . For r 32 > 0.1, the convergence of series in Eq. (j4.2j) is extremely 
good. However, the convergence slows down for the case when r 32 < 0.1. This problem may 
be solved as follows. Applying the identity from Eq. ()3.4jl again, we break the first sum in 
Eq. fl4.2j) in three parts 



G 3d (x 1 ,x 2 ,x 3 ) = G ELC (x 1 ,x 2 ,x 3 ) +G ailh (x 1 ,X2,x 3 ) + { 1 + C j ^ ) J> . ( L-Vi 




where 



9 



Gelc(xi, x 2 , x 3 ) 



1 U 



E 



7T 



exp(-7r7 30 ({m}))cosh 7r7 30 ({m}) (2f^ 



7rhl 2 ^ 73o({^}) 

mi,m2 



sinh [7T73o({m})] 



(4.6) 



x 



J^J cos I 27rmj 



T 



and 



I/3 ^ — \ 7T / / r t \ P^3 1 

Gsiab(xi,x 2 ,a;3) = -77- > — -y-^-exp -27r7 3 o({m})— (4.7) 
xJIcos (^ mi |)-|||* 3 |. 

i=l x ' 

We note an important aspect of this break up of the sum in Eq. ()4.2j) in two parts. Eg. ()4.7j) is 
independent of l 3 as /s/73o({^}) does not depend on / 3 . In fact the expression in Eg. ()4.7j) is 
a three dimensional Coulomb sum for a cell which is open along the x 3 direction and periodic 
along xi and x 2 - Thus the sum in Eg. (J4.7|) corresponds to the slab geometry. Note that the 
suffix ELC stands for the so called "Electrostatic correction term ", a phrase borrowed from 
Ref. |l2j At this point it is worth while to recast the last term in Eg. (j4.5|) in a different 
form, which will prove to be useful later in the discussion. Suppose we have n charges in a 
charge neutral unit cell Yli Qi = 0- Let us assume that the position of the qi is denoted by 
(xu, X21, X3i). Then the third term in Eg. (J4.5|) will give rise to a term in the total energy. 
This term will be given by 



E. 2 « 



l^2qiqj\x 3i -x 3j \ 2 ) , (4.8) 



which after expanding the argument, and using the charge neutrality condition gives 

E z = -yMl (4.9) 

where M 3 = qiX 3i stands for the total dipole moment along the x 3 direction. 

Let us now consider the convergence of Gelc and G s i a b- The function Gelc decays as 
exp (— 27T73o({m})[l — ja^l/^]). Thus we see that Gelc converges exponentially fast for 
< r 3 < 0.5. In fact the slowest convergence of Gelc occurs for the case r 3 = 0.5, but even 
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this slowest convergence goes as exp (— 71730 ({m})), which is extremely fast keeping in mind 
the inequality of Eq. (|4.4|) . 

Now we consider the convergence of G s iab- The previously mentioned problem of conver- 
gence still persists and G s i a b fails to converge fast when < r 32 < 0.1. So, the next step 
is to separate out this diverging behavior towards small value of r 32 . For that purpose we 
break the sum over m/s in Eq. (|4.7|) as follows: 

t 

mi, m 2 m' 2 m [,m 2 

where m' 1 implies that the term corresponding to mi = is not to be included. Thus we 
break up G s iab as 

G B i ab (x 1 ,X2,x 3 ) = Gi(x 2 ,x 3 ) + G 2 (x 1 ,x 2 ,x 3 ), (4.10) 

where 



Gi(x 2 ,x 3 ) = 




(4.11) 



and 



G 2 (xi,x 2 ,x 3 ) 




(4.12) 



First we obtain G\ in a closed form as follows. We may employ the identity from Eq. fl3.6|) 
to obtain 



Gx(x 2 ,x 3 ) 



cosh ( 27T— - 



x 2 

COS I 27T — 

h 



(4.13) 



ln ( 2 ) , ~J**\ 

+ Z7T- 



h hh 

As discussed in the appendix A, G\ has a logarithmic divergence when x 2 /l 2 and x 3 /l 2 tend 
to zero. As we will see soon, a similar logarithmic divergence with opposite sign arises from 
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the term G 2 . These two divergences cancel each other to give a finite contribution to G s i a b 
towards small values of 22 and X3. 

We consider the case of G 2 from Eq. (J4.12|) . Applying the Poisson summation rul©^, the 
sum over m 2 in Eq. (J4.12|) may be transformed to a sum involving Bessel functions of the 
second kincU^: 



1 exp i^-\z\\ja 2 + (x m )^ 



|f|E» ., / exp(2, in »|) (4.14) 

= ^ -K" (a a/z 2 + (2 + 5m)^ 

m 

Identifying 



<5 = Z2, -2 = £3, a = 2vr— - — and x = x 2 , (4-15) 

'1 

we can write 



G 2 (ar 1>a ; 2 ,x 3 ) = j ^ K ^Tr^^fo + m 2 l 2 f + x£\ (4.16) 

17^ ,1712 



x cos (2^^1-^1x31 
The sum in Eq. (j4.16|) may be expressed in two parts as 



G 2 (xi,x 2 ,x 3 ) = j ^ K o \ 2 ' K ~J~\I ( X2 + m ^) 2 + x l J cos ( 27rmi y") ( 4 ' 17 ) 

+ (2^V / ^I + % 2 ) cos (27^!) - ^|x,|. 

We note that the first term in Eq. ()4.17|) has no convergence problem as x 2 and x% are positive 
numbers and l 2 >l\. This term will convergence even for the case when < x 2 and X3 are 
zero. The convergence of G 2 and thus that of G s iab and G3& depend upon the ratio 

^M±£|T, (4 , 8) 
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which appears in the second term on the rhs of Eq. (j4.17|) . For p > 0.1, Eq. (j4.17|) will have 
a very good convergence. However, if x 2 and 23 are such that the condition p > 0.1 is not 
satisfied then we should transform Eq. (j4.17|) further. This can be done by using the results 
derived in appendix B where it is shown that 




(4.19) 



2 7 {i/j(N + x) + ip(N -x)} 
~h h 



+ fE( 1 / 2 )p 2 '(C(2/ + l,iV + x)+C(2/ + l,iV-x)), 
1 1=1 ^ ' 

where x = x\jl\ and ip and ( stand for digamma and Hurwitz Zeta function respectively. 
iV > 1 is the smallest integer satisfying the condition N > p + x. Thus we can choose N = 1, 
as even for the worst case one has p — 0.1 and x = 0.5 . However, for better convergence it 
is desirable that one chooses N such that N > p + 1. 

We can now write the following short algorithm to calculate G s i a b • First we set our axis 
such that l 3 > l 2 > l\. Next, using the periodic boundary conditions, the separation between 
two particles can always be reduced in such a way that the individual components satify 
< Xi < li. Thus, the values of rj = Xi/li lie between and 1. From the inherent symmetry 
of the problem, energy corresponding to eight different separations of (^4p, , is 
the same. This essentially means that we can concentrate our attention on only on those 
separations between the particles which correspond to < r,i < 0.5. If some > 0.5, we can 
replace it with 1 — ?y Next, we look at the value of = ^z/h- If ^32 > 0.1, we can combine 
Eq. (j4.13ft with Eq. (j4.12|) to obtain the following form for G s i a b 
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(T ln 



X 3 



u 



cosh 2ix— — cos 2n— 



x 2 



(4.20) 



ln(2) + 1 1 



l\ 7T Z1Z2 



m, ,m2 



71 



mi \ 1 Z m2 
(1 / V ' 2 



-W(¥)*(¥HH-(*-f) 



However if < r 32 < 0.1, then we look at the value of p, which is defined in Eq. (|4.18|) . If 
p > 0.1, we should use the following form of G s iab which is obtained after combining Eqs 
flUIH and (14371) : 



G a l ah (xx,X 2} X 3 ) = - J- l n 
'l 



^3 



cosh 2tt— — cos 27T— 



X2 



(4.21) 



^ + \ E A ° (2vr^^ 2 + m 2 / 2 ) 2 + x|) cos ^vrm^ 



mi=l 



/ 2ivm\ 



(x 2 + X3J 



2\V2 



COS 



/ 27rm! 



If p < 0.1 then we use the identity in Eq. (j4.19|) to write G s i a b as 



G , s iab(xi,x 2 ,x 3 ) = -—In 



X3 \ / %2 

cosh I 27T— I — COS I 27T— 



(4.22) 



-^ + ^^» (2^^ + m^ + xl) cos (2*^!) 



+ - <( 7 + In 



(^ + x|) 1/2 ~ 



JV-l 



2/1 



1 



+ 



•y^ ~\~ x\ -\- X3 



1 



rti=l 



\j p 2 + (ni + x) 2 \J p 2 + ( n i - ~%y 

2 7 {ip(N + x) + ip(N-x)} 



h 



p 21 (C (2l + l,N + x) + £ (21 + 1,N — x)) 



Although Eq. (|4.22|) is meant to be used only when p < 0.1, the equation is defined for all 
values of p as long as N is chosen such that N > p + 1. The series given in Eq. (j4.22|) is valid 
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when both %2 and £3 are non zero. In this case, the argument of the first logarithmic term 
on the rhs of Eq. (|4.22j) is always greater than zero. However, for very small values of £2 and 
£3 (say both less than e = 10 -3 ) the first and the fourth terms diverge. In such situation 
one should combine the diverging terms together using the function L defined in appendix 
B. 

We have thus shown how to compute G s iab for all regions of the unit cell. Similar results 
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15 and 16. Results 



for the slab geometry have previously been obtained in Refs. 
in Eqs. (j4.2U|) and (J4.22JI correspond respectively to "near" and "far" formulae derived by 
Arnold et a/.—. Also, it is an easy matter now to obtain expressions for G^d from Eq. (j4.5j) . 
One can obtain the self energy for a 3D system as, 



G™ { = ]im I G M ( Xl ,X2,x 3 ) --====) (4.23) 

(»i,x 2 ,a;3)->(0,0,0) V X \ + x\ + x\ 



\h_ 77 exp (-7r7 30 ({m})) 

n ll12 mT^ 2 73o({m}) sinh(7T73o({m})) 



. JjA h 7T 2 fAnlA , 2 7 



nij ,m 2 

where for G 3d we use Eqs. (J4.5j) and (|4.22jl . 

V. MADELUNG CONSTANTS 

Using the formulas developed above, it is an easy matter to obtain expressions for the 
Madelung constants of NaCl and CsCl. A simple structural analysis of CsCl easily leads to 
expression, 



^ A v< 1 - G 3d (-,-,-)- G 3 ^, ( 5. 1 ) 



and similarly for NaCl we see 



^NaCl — ^ 



G:u I ~ I, I) + 3G3H (V 0, 1} - 3G sri I ~ i ) - G 3d 



3d ' 2' 2' 2 7 d V 2 7 ~ V 2 2 ' ~ ! 11 



(5.2) 

where Eq. f)4.5j) can be used for G^{x\, X2, £3) with all k's set equal to one. From the above 



equations we obtain the following expressions for Madelung constants of CsCl and NaCl: 
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Mi 



CsCl 



7T 



exp l —it 



7T Z — ' -v/m? + 
mi, ma V 1 1 



mi 



sinh ( 7r y mf + mf) 



(5.3) 



iV (27r |m 1 m 2 |) — In (47r) + 7 — it 



and 



1 ' vr (exp \-it^m\ + m|l - (-1)™^ - 3 + 3 (-1 

2 M NaC1 = V ' J 

7T ^— ' 



mi 



71 mi ,m 2 V^i + m 2 sinh(7r v /m|Tm|) 

(5.4) 

2 I ^ K (27r |mim 2 |) - In (Ait) + 7 - vr 

V™'l> m 2 

Restricting the sum over mi and m 2 between —4 and +4, a simple calculation on Mathemat- 
ica gives a M CsC1 value correct up to 1CT 8 and M NaC1 value correct up to 1CT 6 . In addition 
we also obtain a simple relationship between the two Madelung constants: 



esc I Try (2mi + l) 2 + m 2 
2 M NaC i = M CsC i + 6 / = 7 • ( 5 - 5 ) 

mi,m 2 y (2m 1 + 1) + m| 



This interesting relationship was first established by Hautolj^ in seventies using Hankel 
integrals and Schloimilch series. 



VI. CONCLUSION 



Complete expressions for Coulomb sum for a rectangular cell in 2D and an orthorhombic 
cell in 3D were derived. We also obtained expressions for the self energies. The expressions 
obtained provide convergence in all parts of the unit cell. Considerable simplification has 
been achieved over Sperb's work 2 - in terms of deriving the equations. The proposed formula 
for the potential energy when the two charges are very close, differs from that of Sperb. In 
particular, when the charges are close together, Sperb's 2 formula has a triple sum (Eq. 2.4 
and 2.7). In our expression, we have at most a double sum. Similar results for 3D case have 
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previously been obtained by Strebel using a rather involved procedure^. Our results do not 
require any convergence parameter like that used in Ewald sums, neither do our formulas 
involve any complementary error functions. These error functions in an Ewald sum are a 
source of loss of precision when calculating Madelung constants to higher accuracies. 

In retrospect, we see that these results could be derived in another way by starting off 
with the Green function expression for the 2D + h slab geometry system and then adding 
the ELC term which takes into account the rest of the layers. This way we will get only the 
first two terms of Eq. (j4.5j) . The third term is then obtained by adding a term proportional 
to M| from outside, where M3 stands for the component of the total dipole moment along 
the X3 direction. In the present work this dipole term arises naturally, as shown in Eq. 
(14. 9|) . This dipole term has been discussed by Smith^. Thus this slabwise summation plus 
a term dipole term added from outside, apart from an unimportant constant, leads to the 
same expression as in Eq. (j4.5l) . Thus, our Eqs. (j4.5|) and ()4.9|) can be viewed as an alternative 
derivation of Eq.(4) in Ref. Il2|. 

An advantage of the method developed here is that one can achieve better time scaling 
in a simulation. Using the expressions presented in this paper, the time to calculate forces 
and energy for a 3D system in a computer simulation scales as N 2 , where N is the number 
of charges present in the unit cell. However, one can achieve N 5 ^ 3 ln(iV) 2 scaling after a 
little modification in the expressions presented here. This is the same scaling as achieved by 
Arnold et a/.— for 2D + h system. The scaling remains the same for the two cases because 
the electrostatic correction term can be computed linearly if we remove the contribution of 
the first two closest layers enclosing the unit cell in a given direction as opposed to removing 
the contribution of just one layer as done by Arnold et a/.— and in this paper. Also the 
results presented here can be generalized to a rhombic cell in 2D and a triclinic cell in 3D 20 

Our proposed expressions can be applied to calculation of Madelung constants in 3D. 
Results obtained for the Madelung constants of CsCl and NaCl match with those in the 
literature. 

In conclusion we have provided a very simple derivation of complicated results previously 
obtained by many authors using different, sometimes complicated, techniques ^ 4 ' 13 > 15 i 16 . 
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APPENDIX A: LOGARITHMIC DIVERGENCE 



Consider the function 



L(x, y) = In [cosh?/ — cos re] — In 



y 2 + x 2 



(Al) 



We want to examine the limiting value of L as x and y tend to zero. For this reason we 
expand the argument of the first logarithmic term in Eq. (|A1J) . 

V + x e 



cosh y — cos x 



+ 



y 2 + x 2 
~~ 2! 

y 8 — x 8 



4 4 

y — x 
4! 



6! 



(A2) 



+ o[* 1( V ]. 



Factoring out the first term on the right hand side, Eq. ()A2|1 can be written as 
coshy — cosx = 



9 9 

y z + ar 



2! 



l + |(y 2 -x 2 ) + |(y 4 -a;V + x 4 ) 



(A3) 



+ |(*/ 4 + * 4 ) (y 2 -x 2 )+0[x 8 ,y*]y 



Thus L can be written as 



L(x, y) = In 1 1 + | {y 2 - x 2 ) + | (y 4 - x 2 y 2 + x 4 ) 
+| (l/ 4 + ^ 4 ) (y 2 -x 2 )+0[x 8 !2 / 8 ]}. 



(A4) 



Using the results from Eqs. (jAl|) and (jA4|) in Eq. (j4.13J) . we see that for small values of X2/I2 
and X3//2, Gi can be written as 



Gi(x 2 ,x 3 ) 



In 



2tt' 



, (x 2 + x|) 



/ 2 



ln ( 2 ) , o^l^l 
~ ^ — + 27r TT 



(A5) 



_ — T ( -L ^1 



which clearly shows a logarithmic divergence as X2/I2 and x^/k tend to zero. 
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1. Line Charge 

We commence with the identity^, 



where 



„, .4 /27rmi , 2 2\i/2\ f 2irrrii \ 

f {x 1 ,x 2 ,x 3 ) = — K I — - — {x 2 + x 3 ) Jcosl — - — xi I (A6) 

^ mi=l ^ ^ / V 1 / 

= JT I 7 + ln ) | + TWTWT^ + s (I1 ' I2 ' I3) • 



S (xi, x 2 , x 3 ) = \ i 2 ( A7 ) 

\jx\ + x§ + (nZi - xi) 2 



n=l 



+ , — - - • (A8) 

ya;| + x| + (nZi +xi) 2 n / 

We can further transform the identity in Eq. ()A6j) along the lines worked out by Strebel^ 
and Arnold et al^. Let us look at 



oo 



n=N 



1 



oo 



h 



n=l 



l 



2 n 



\ 



2 n + N - 1 



1 oo 



1 Wrt 1 - 



n=i \\l p 2 + (n + y) J n=l 
where N > 1, y = iV — 1 + xi/l\ and 

- o . oil/2 



p " /l ' h 



(A9) 



(A10) 



(All) 



(A12) 



Assuming p < |1 + y\, the Binomial expansion of the first term in the above equation gives 
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Y(~ 1/2 )p 2p — o— TT (A13) 

P J \n + y\ 2p+1 V ' 



sjp* + (n + yf P ' \ n + y\ 

h \ * ) \n + y\ 2p+1 + 



n + y\' 



where ( stands for the Binomial coefficient. We can take the sum over n inside and 
obtain 



^-ixi~Ty%]^ <ai4) 

p=l x ' n=l I <* I 

+ r> 1 ,-- +r) -• A15 

n=l x 7 n=l 



Now, using the definition of the Hurwitz Zeta function, 

oo 1 

fc =o ( fc + 

we obtain 

E l 7 |2p+1 =C(2p+l,l + y). (A17) 

Also the second sum in Eq. (|A14|) is easy to obtain. By the definition of the digamma function 
ip we have 

Efi-4-r--) =-7-^(i + i/). (A18) 

Thus /i (p, Xi) can be written as 

k(ftIl) = ^_^ + if;(-^) A(a + 1 , 1 + B) + ^I (A19) 

1 1 1 1=1 v 7 1 n=l 

Using Eqs. (lAll) . (jXTTjt and (lAlSl) we obtain, 
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1 JV " 1 / 1 1 2 \ 

S ( Xl , x 2 , xs) = 77 E| / + / - - I ( A2 °) 

\J p 2 + (n - xi) 




2 



\J p 2 + (n- xi) 



2 n 



1 



'p 2 + (n + xi) Jp 2 + (n — Xi j 
27 ip(N + xi) + ip(N-xi) 



a/x? + xj + X 



+ r£( 1 / 2 )p 2 '(C(2/ + l,iV + x) + C(2/ + l,iV-a;)), 
1 i=i ^ ' 

Note that for Eq. (|A20J) to be valid, the condition is that p < \l+y\, where y = N— l±xi//i. 
Keeping in mind that Xi > we get p < \N ± Xi/Zi|, which will be satisfied if N > p + x. 
Combining Eq. (|A6|) and Eq. ()A20|) gives us the result that we set out to prove. 
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